* Figure: RANDOMIZATION INFERENCE

use "${output}panel_r0_r1_r2.dta", clear // Load merged analysis panel

// Set seed
set seed 16523

********************************************************************************

* Keep households surveyed in each round
keep if sample_household == 1

// Reposition baseline value of household controls
foreach var of varlist household_size children_under_five {
	gen b_`var' = `var' if surveyround == 0
	bys hh_id (surveyround) : replace b_`var' = b_`var'[1] if mi(b_`var')
	}

********************************************************************************

// Analyses
keep if surveyround == 1

// "True" value
qui reghdfe purchased_intervention_stove i.treatment##i.chirag_strata b_*, absorb(districtcode) vce(cluster uniquegrp)
local true_effect = _b[1.treatment#1.chirag_strata]
local true_effect : di %5.2f `true_effect'
disp "`true_effect'"

// Program for random assignment
capture program drop ngo_inference
program define ngo_inference, rclass
	preserve	
		sort hh_id
		bys gp_code (hh_id) : gen rand = runiform() if _n == 1
		gen chirag_strata_random = .
		replace chirag_strata_random = 1 if rand > 0.5 & !mi(rand)
		bys gp_code (chirag_strata_random) : replace chirag_strata_random = chirag_strata_random[1] if mi(chirag_strata_random)
		replace chirag_strata_random = 0 if mi(chirag_strata_random)		
		reghdfe purchased_intervention_stove i.treatment##i.chirag_strata_random b_*, absorb(districtcode) vce(cluster uniquegrp)
		return scalar beta_ngo_inference = _b[1.treatment#1.chirag_strata_random]
	restore
end

// Simulation
qui simulate beta_ngo_inference_results = r(beta_ngo_inference), reps(1000) : ngo_inference

// Percentiles
_pctile beta_ngo_inference_results, nq(1000)
local p25 = `r(r25)'
local p5 = `r(r50)'
local p950 = `r(r950)'
local p975 = `r(r975)'

// Graph
twoway (kdensity beta_ngo_inference_results, range(`p25' `p975') color(gs14) recast(area)) ///
	|| (kdensity beta_ngo_inference_results, range(`p5' `p950') color(gs11) recast(area)) ///
	|| (kdensity beta_ngo_inference_results, ///
		lstyle(solid) ///
		ytitle(Density) ///
		xtitle("{&beta}(TREATMENT{sub:j} × NGO{sub:j})")) ///
	|| (scatteri 0 `true_effect' 5.8 `true_effect' , recast(line) lpattern(dash)) ///
	|| (scatteri 5.8 `true_effect' "Estimate = `true_effect'", msymbol(none) mlabsize(*1.25) mlabposition(12) ///
			legend(order(2 1) label(1 "95% C.I.") label(2 "90% C.I.")))

graph export "${results}figure_ngo_randomization_inference.png", ///
	width(5000) replace

gr drop _all

// Check for share of large results
gen ngo_inference_p = (beta_ngo_inference_results > `true_effect')
sum ngo_inference_p
